@eclare108213 : For the snwredist options, it seems we're only supporting ITDrdg and none. Are bulk, 30percent and ITDsd defunct?

We're supporting ITDrdg, bulk and none. The old 30percent is available - it's set using 'bulk' with a (new) namelist parameter that can be any percentage. ITDsd is defunct.

Good to know. I'll have to look into how the bulk option would work in mpas-seaice. Does 'bulk' use the effective snow density?

I think bulk works the same way as 30percent did before, just generalized the percentage. So it only rearranges the snow without changing the density. And yes the density tracer is still available in that case - it should remain constant, although I haven't looked at that yet after this latest round of changes.

I'm concerned that the way it's implemented, we are requiring the density tracer for the bulk option. Otherwise it will fail. And it doesn't do anything.

It doesn't fail in Icepack (or CICE). Can MPAS-SI not implement it that way? With it, analysis scripts looking for it (for comparison with ITDrdg) won't fail, and its values can be verified to be constant.
Have you tried running this icepack version in your E3SM tests? I think it's coded so that if the driver doesn't have the tracer, icepack will still handle it gracefully. But I'm not sure about that. If it doesn't work, then I can change it.

I haven't tested (recently) snwredist = bulk and snow effective density off, but I will. But tracers are passed to icepack in the tracer array which in mpas is only made up of active tracers. The index for nonactive tracers is set to 0. So in something like

do k = nt_rhos, nt_rhos+nslyr-1
trcrn(k,n) = max(trcrn(k,n)-rhosmin, c0)

I suspect F90 won't like that nt_rhos = 0, and the rest of the code will overwrite something valuable.

I could be wrong. I've almost got your changes implemented. I'll try running with different options.

The changes 141, 267 and 268 would be a problem with other ktherm options and the snow grain tracer active. Since the internal assignments of massice and massliq are wrapped in a ktherm == 2 statement, they would be sent as zeros to thickness_changes.

My first look didn't catch this. Probably the simplest fix is to just pull the assignment statements out of the ktherm == 2 if-else and move to after both the temperature_changes routines.

Do we need to support advanced snow physics on BL99 thermo? We could abort if that combination is attempted. Or we could initialize massice and massliq according to the BL99 assumptions (use rhos for massice and 0 for massliq?).

I think just moving it out of the if statement is easiest. It has it's own if snwgrain ... wrapper. I'm certainly not in favor of supporting a massive number of config combinations, but the snow stuff doesn't really care about mushy or BL99. That's not to say I've tested it with BL99...

Do you mean line 333ff?

Yes

Okay, I agree that should be changed.

@eclare108213 : line 1294 in icepack_tracers.F90 should also be changed to rhosnew.

done

This else part of this statement in icepack_snow.F90 is a problem. rhos_cmpn can't be used at all if snow effective density is false otherwise it will overwrite ice area.
Lines 318 to 320
else
rhos_cmpn(:,:) = rhos
endif

And the first part of the if statement
if (trim(snwredist) /= 'none') then
should be changed to
if(snwredist(1:3) = 'ITD') then

Agreed. Will it work to just remove the else part? This is part of the reason for defining the all of the tracers when tr_snow=T.

yes, deleting L318 and 319 works

@eclare108213 : This expression (L1245 and 1246) forces the ice mass tracer to maintain it's current ice density by adding/subtracting mass from the snow liquid tracer.

1240  if (hsn > puny) then    ! add snow with enthalpy zqsn(1)                                                            
1241        dhs = econ / (zqsn(1) - rhos*Lvap) ! econ < 0, dhs > 0  
1242         mass  = massice(1) + massliq(1)                                                                                  
1243         massi = c0                                                                                                       
1244         if (dzs(1) > puny) massi = c1 + dhs/dzs(1)
1245         massice(1) = massice(1) * massi                                                                                  
1246        massliq(1) = max(c0, mass + rhos*dhs - massice(1)) ! conserve new total mass 

It's only going to be round off level changes at this point because the ice mass is always near rho, but it's not what we want. It should just be

massice(1) = massice(1) + dhs * rhos

[edited formatting only]

To clarify, @njeffery, you are suggesting that we should only change massice(1) and not change massliq(1)? Having thought through this, I think you are right. It's pretty confusing, though, so I'm recording my thought process here for later reference:
Condensation could be done either way, with all of the condensation becoming solid ice or dividing it between ice and liquid. I think the massi code above divides it evenly based on what the tracers were to begin with. But the latent heat values would probably be different in that case. Is that the issue here? This calculation is using Lvaporization, which is between liquid and gas, but if we're condensing from gas directly to ice, shouldn't it be Lsublimation? Answer: No, because the enthalpy definition that we use in the model includes the heat needed to convert between solid and liquid phases (and to bring the liquid to 0 C). Here, we are only changing the ice and liquid mass tracers, which shouldn't affect the solution at this point in the code -- the liquid mass may be used later for melt ponds, which would change the solution in those configurations.

@eclare108213 : yes. Just massice. New snow formation has massice = rhos and this is most consistent with both our energy and mass exchanges. The points you're making above about vapor to liquid and vapor to ice are also worth thinking more about, especially if we ever start changing snow density in the thermodynamics. But we okay for conservation. You'll see differences through the melt ponds and the snow grain radius.